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ABSTRACT 

We present the angular power spectrum derived from the first-year Wilkinson Microwave Anisotropy 
Probe (WMAP) sky maps. We study a variety of power spectrum estimation methods and data combina- 
tions and demonstrate that the results are robust. The data are modestly contaminated by diffuse Galactic 
foreground emission, but we show that a simple Galactic template model is sufficient to remove the signal. 
Point sources produce a modest contamination in the low frequency data. After masking ~700 known 
bright sources from the maps, we estimate residual sources contribute ~3500 fiK 2 at 41 GHz, and ~130 
pK 2 at 94 GHz, to the power spectrum [/(/+ 1)C//27 t] at l = 1000. Systematic errors are negligible com- 
pared to the (modest) level of foreground emission . Our best estimate of the power spectrum is derived 
from 28 cross-power spectra of statistically independent channels. The final spectrum is essentially inde- 
pendent of the noise properties of an individual radiometer. The resulting spectrum provides a definitive 
measurement of the CMB power spectrum, with uncertainties limited by cosmic variance, up to / ^ 350. 
The spectrum clearly exhibits a first acoustic peak at l = 220 and a second acoustic peak at / ~ 540 (Page 
et al. 2003b), and it provides strong support for adiabatic initial conditions (Spergel et al. 2003). Kogut 
et al. (2003) analyze the Cj E power spectrum, and present evidence for a relatively high optical depth, 
and an early period of cosmic reionization. Among other things, this implies that the temperature power 
spectrum has been suppressed by ~30% on degree angular scales, due to secondary scattering. 
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1. INTRODUCTION 

The Wilkinson Microwave Anisotropy Probe (WMAP) mission was designed to measure the CMB anisotropy 
with unprecedented precision and accuracy on angular scales from the full sky to several arc minutes by producing 
maps at five frequencies from 23 to 94 GHz. The WMAP satellite mission (Bennett et al. 2003a) employs a matched 
pair of 1.4m telescopes (Page et al. 2003c) observing two areas on the sky separated by ~ 141°. A differential ra- 
diometer (Jarosik et al. 2003a) with a total of 10 feeds for each of the two sets of optics (Barnes et al. 2003; Page 
et al. 2003c) measures the difference in sky brightness between the two sky pixels. The satellite is deployed at the 
Earth-Sun Lagrange point, Z> 2 , and observes the sky with a compound spin and precession that covers the full sky 
every six months. The differential data are processed on the ground to produce full sky maps of the CMB anisotropy 
(Hinshaw et al. 2003b). 

Full sky maps provide the smallest record of the CMB anisotropy without loss of information. They permit a 
wide variety of statistics to be computed from the data - one of the most fundamental is the angular power spectrum of 
the CMB. Indeed, if the temperature fluctuations are Gaussian, with random phase, then the angular power spectrum 
provides a complete description of the statistical properties of the CMB. Komatsu et al. (2003) have analyzed the first- 
year WMAP sky maps to search for evidence of non-Gaussianity and find none, aside from a modest level of point 
source contamination which we account for in this paper. Thus, the measured power spectrum may be compared to 
predictions of cosmological models to develop constraints on model parameters. 

This paper presents the angular power spectrum obtained from the first-year WMAP sky maps. Companion 
papers present the maps and an overview of the basic results (Bennett et al. 2003b), and describe the foreground 
removal process that precedes the power spectrum analysis (Bennett et al. 2003c). Spergel et al. (2003), Verde et al. 
(2003), Peiris et al. (2003), Kogut et al. (2003), and Page et al. (2003b) discuss the implications of the WMAP power 
spectrum for cosmological parameters and carry out a joint analysis of the WMAP spectrum together with other CMB 
data and data from large-scale structure probes. Hinshaw et al. (2003b), Jarosik et al. (2003b), Page et al. (2003a), 
Barnes et al. (2003), and Limon et al. (2003) discuss the data processing, the radiometer performance, the instrument 
beam characteristics and the spacecraft in-orbit performance, respectively. 

A sky map AT(n) defined over the full sky can be decomposed in spherical harmonics 

/ 

A7’(n) = ^^a ;m y, m (n) (1) 

/>0 m=—l 

with 

a lm = J dn n AT(n)Y; m (n), (2) 

where n is a unit direction vector, and T/ m (n) are the spherical harmonic functions evaluated in the direction n. 

If the CMB temperature fluctuation AT is Gaussian distributed, then each is an independent Gaussian deviate 

with 

<«/«)= 0, (3) 

and 

) ~ $iv <W Q , (4) 

where Q is the ensemble average power spectrum predicted by models, and <5 is the Kronecker symbol. The actual 
power spectrum realized in our sky is 

, m=-l 


( 5 ) 
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In the absence of noise, and with full sky coverage, the right hand side of equation (5) provides an unbiased estimate 
of the underlying theoretical power spectrum, which is limited only by cosmic variance. However, realistic CMB 
anisotropy measurements contain noise and other sources of error that cause the quadratic estimator in equation (5) to 
be biased. In addition, while WMAP measures the anisotropy over the full sky, the data near the Galactic plane are 
sufficiently contaminated by foreground emission that only a portion of the sky (~85%) can be used for CMB power 
spectrum estimation. Thus, the integral in equation (2) cannot be evaluated as such, and other methods must be found 
to estimate Q. 

In Appendix A, we review two methods that have appeared in the literature for estimating the angular power 
spectrum in the presence of instrument noise and sky cuts. The first (Hi von et al. 2002) is a quadratic estimator that 
evaluates equation (2) on the cut sky yielding a “pseudo power spectrum” C/. The ensemble average of this quantity is 
related to the true power spectrum, Q by means of a mode coupling matrix Gw (Hauser & Peebles 1973). The second 
method (Oh et al. 1999) uses a maximum likelihood approach optimized for fast evaluation with WMAP-like data. In 
Appendix A we demonstrate that the two methods produce consistent results. 

A quadratic estimator offers the possibility of computing both an “auto-power” spectrum, proportional to |a/ m | 2 , 

and a “cross-power” spectrum, proportional to where the a\ m coefficients are estimated from two indepen- 

dent CMB maps, i and j. This latter form has the advantage that, if the noise in the two maps is uncorrelated, the 
quadratic estimator is not biased by the noise. For all cosmological analyses, we use only the cross-power spectra be- 
tween statistically independent channels. As a result, the angular spectra are, for all intents and purposes, independent 
of the noise properties of an individual radiometer. This is analogous to interferometric data which exhibits a high 
degree of immunity to systematic errors. The precise form of the estimator we use is given in Appendix A. 

The plan of this paper is as follows. In §2 we Teview the properties of the WMAP instrument and how they 
affect the derived power spectrum. In §3 we present results for the angular power spectra obtained from individual 
pairs of radiometers, the cross-power spectra, and examine numerous consistency tests. In §4, in preparation for 
generating a final combined power spectrum, we present the full covariance matrix of the cross-power spectra. In §5 
we present the methodology used to produce the final combined power spectrum and its covariance matrix. In §6 we 
compare the WMAP first-year power spectrum to a compilation of previous CMB measurements and to a prediction 
based on a combination of previous CMB data and the 2dFGRS data. We summarize our results in §7 and outline 
the power spectrum data products being made available through the Legacy Archive for Microwave Background Data 
Analysis (LAMBDA). Appendix A reviews two methods used to estimate the angular power spectrum from CMB 
maps. Appendix B describes how we account for point source contamination. Appendix C presents our approach to 
combining multi-channel data. Appendix D describes how the foreground mask correlates multipole moments in the 
Fisher matrix, and Appendix E collects some useful properties of the spherical harmonics. 


2. INSTRUMENTAL PROPERTIES 

The WMAP instrument is composed of 10 “differencing assemblies” (DAs) spanning 5 frequencies from 23 
to 94 GHz (Bennett et al. 2003a). The 2 lowest frequency bands (K and Ka) are primarily Galactic foreground 
monitors, while the 3 highest (Q, V, and W) are primarily cosmological bands. There are 8 high frequency differencing 
assemblies: Ql, Q2, VI, V2, and W1 through W4. Each DA is formed from two differential radiometers which are 
sensitive to orthogonal linear polarization modes; the radiometers are designated 1 or 2 (e.g., VI 1 or W12) depending 
on which polarization mode is being sensed. 

The temperature measured on the sky is modified by the properties of the instrument. The most important prop- 
erties that affect the angular power spectrum are finite resolution and instrument noise. Let C\ l denote the auto or 
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cross-power spectrum evaluated from two sky maps, i and i f , where i is a DA index. Further, define the shorthand 
i = (i,z') to denote a pair of indices, e.g., (Q1,V2). This spectrum will have the form 

C^w^+Ni, (6) 

where w) = b\ b\ pj is the window function that describes the combined smoothing effects of the beam and the finite 
sky map pixel size. Here b\ is the beam transfer function for DA /, given by Page et al. (2003a) [note that they reserve 
the term “beam window function” for (b|) 2 ], and pi is the pixel transfer function supplied with the HEALPix package 
(Gorski et al. 1998). N j is the noise spectrum realized in this particular measurement. On average, the observed 
spectrum estimates the underlying power spectrum, Q, 

(C})=w{C/ + (A/)fe, (7) 

where (N\) is the average noise power spectrum for differencing assembly i, and the Kronecker symbol indicates that 
the noise is uncorrelated between differencing assemblies. To estimate the underlying power spectrum on the sky, C/, 
the effects of the noise bias and beam convolution must be removed. The determination of transfer functions and noise 
properties are thus critical components of any CMB experiment. 

In §2.1 we summarize the results of Page et al. (2003a) on the WMAP window functions and their uncertainties. 
We propagate these uncertainties through to the final Fisher matrix for the angular power spectrum. In §2.2 we present 
a model of the WMAP noise properties appropriate to power spectrum evaluation. For cross power spectra (i / i l 
above), the noise bias term drops out of equation (7) if the noise between the two DAs is uncorrelated. These cross- 
power spectra provide a nearly optimal estimate of the true power spectrum, essentially independent of errors in the 
noise model, thus we use them exclusively in our final power spectrum estimate. The noise model is primarily used to 
propagate noise errors through the analysis, and to test a variety of different power spectrum estimates for consistency 
with the combined cross-power spectrum. 


2.1. Window Functions 

As discussed in Page et al. (2003a), the instrument beam response was mapped in flight using observations of 
the planet Jupiter. The signal to noise ratio is such that the response, relative to the peak of the beam, is measured to 
approximately -35 dB in W band, the band with the highest angular resolution. The beam widths, measured in flight, 
range from 0?82 at K band down to 0?20 in some of the W band channels (FWHM). Maps of the full two-dimensional 
beam response are presented in Page et al. (2003a), and are available with the WMAP first-year data release. The 
radial beam profiles obtained from these maps have been fit to a model consisting of a sum of Hermite polynomials 
that accurately characterize the main Gaussian lobe and small deviations from it. The model profiles are then Legendre 
transformed to obtain the beam transfer functions b\ for each DA i. Full details of this procedure are presented in Page 
et al. (2003a), and the resulting transfer functions are also provided in the first-year data release. We have chosen to 
normalize the transfer function to 1 at / = 1 because WMAP calibrates its intensity response using the modulation of 
the CMB dipole (/ = 1). This effectively partitions calibration uncertainty from window function uncertainties. 

The beam processing described above provides a straightforward means of propagating the noise uncertainty 
directly from the time-ordered data through to the final transfer functions. The result is the covariance matrix Ej, w 
for the normalized transfer function. Plots of the diagonal elements of E‘ bir are presented in Page et al. (2003a). 
The fractional error in the transfer functions b\ are typically 1-2% in amplitude. In the end, these window function 
uncertainties dominate the small off-diagonal elements of the final covariance matrix for the combined power spectrum 
(see §4). Additional observations of Jupiter will reduce these uncertainties. 
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An additional source of error in our treatment of the beam response arises from non-circularity of the main beam. 
The effects of this non-circularity are mitigated by WMAF s scan strategy which results in most sky pixels being 
observed over a wide range of azimuth angles. The effective beam response on the sky is thus largely symmetrized. We 
estimate that the effects of imperfect symmetrization produce window function errors of < 1 % relative to a perfectly 
symmetrized beam window function (Page et al. 2003a; Hinshaw et al. 2003b). This error is well within the formal 
uncertainty given in Kw In §5 we infer the optimal power spectrum by combining the 28 cross-power spectra 
measured by the 8 high frequency DAs Q1 through W4. As part of this process we marginalize over the window 
function uncertainty, which automatically propagates these errors into the final covariance matrix for the combined 
power spectrum. Both the combined power spectrum and the corresponding Fisher matrix are part of the first-year 
data release. 


2.2. Instrument Noise Properties 

The noise bias term in equation (7) is the noise per ai m mode on the sky. If auto-power spectra are used in the final 
power spectrum estimate, the noise bias term must be known very accurately because it exponentially, dominates the 
convolved power spectrum at high /. If only cross correlations are used, the noise bias is only required for propagating 
errors. Our final best spectrum is based only on cross correlations, and is independent of this term. However, as an 
independent check of our results, we evaluate the maximum likelihood spectrum based on a combined Q+V+W sky 
map. The noise bias must be estimated accurately for this application. 

In the limit that the time-ordered instrument noise is white, the noise bias is a constant, independent of l. If the 
time-ordered noise has a 1 // component, the bias term will rise at low l. In this subsection we estimate the noise bias 
properties for each of the high frequency WMAP radiometers based on the time-ordered noise properties presented in 
Jarosik et al. (2003b). While the WMAP radiometer noise is nearly white by design (Jarosik et al. 2003a) with 1 // 
knee frequencies of less than 1 0 mHz for 9 out of 1 0 differencing assemblies, one of the radiometers (W4 1 ) has a 1 // 
knee frequency of ~45 mHz. The latter is large enough that the deviations of (Nj) from a constant must be accounted 
for. 


The most reliable way to estimate the effects of 1 // noise on the measured power spectra is by Monte Carlo 
simulation. Using the pipeline simulator discussed in Hinshaw et al. (2003b) we have generated a library of noise maps 
with flight-like properties. Specifically we have included flight-like 1 / / noise in the simulated time-ordered data, and 
have run each full-year realization through the map-making pipeline, including the baseline pre- whitening discussed 
in Hinshaw et al. (2003b). We evaluate the power spectra of these maps using the quadratic estimator described in 
Appendix A with 3 different pixel weighting schemes. (See Appendix A. 1.2 for definitions of the weights, and the l 
range in which each is used.) We define the effective noise as a function of l based on fits to these Monte Carlo noise 
spectra. For the analyses in this paper, we fit the spectra to a model of the form 


In (N}) = 


fimax 

£4 (In/)" 


L n=0 


-1 


( 8 ) 


where the c l n are fit coefficients given in Table 1, with n ma x = 3 for / < 200, and n ^ = 1 for / > 200. 

Figure 1 shows the noise spectrum derived from the simulations for each of the 8 high frequency DAs, using 
uniform weighting over the entire l range. For comparison, we also plot an estimate of the CMB power spectrum from 
§5 in grey. Note that the W4 spectrum is the only one of this set to exhibit deviations from white noise in an l range 
where the signal-to-noise is relatively low, and we believe this simulation slightly over-estimates the 1 // noise in the 
flight W4 differencing assembly (Hinshaw et al. 2003a). 
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2.3. Systematic Errors 

Hinshaw et al. (2003a) present limits on systematic errors in the first-year sky maps. They consider the effects of 
absolute and relative calibration errors, artifacts induced by environmental disturbances (thermal and electrical), errors 
from the map-making process, pointing errors, and other miscellaneous effects. The combined errors due to relative 
calibration errors, environmental effects, and map-making errors are limited to < 15 jjK 2 (2a) in the quadrupole 
moment C 2 in any of the 8 high-frequency DAs. Tighter limits are placed on higher-order moments. We conservatively 
estimate the absolute calibration uncertainty in the first-year WMAPdata to be 0.5%. 

Random pointing errors are accounted for in the beam mapping procedure; the beam transfer functions presented 
by Page et al. (2003a) incorporate random pointing errors automatically. A systematic pointing error of ~10" at the 
spin period is suspected in the quaternion solution that defines the spacecraft pointing. This is much smaller than 
the smallest beam width (~12' at W band), and we estimate that it would produce <1% error in the angular power 
spectrum at l = 1000, thus we do not attempt to correct for this effect. Barnes et al. (2003) place limits on spurious 
contributions due to stray light pickup through the far sidelobes of the instrument. They place limits of <10 p K 2 on 
spurious contributions to Q, at Q through W band, due to far sidelobe pickup. 

A detailed model of Galactic foreground emission based on the first-year WMAP data is presented by Bennett 
et al. (2003c) and is summarized in §3.1. We show that diffuse foreground emission is a modest source of contamina- 
tion at large angular scales (22°). Systematic errors on these angular scales are negligible compared to the (modest) 
level of foreground emission. On smaller angular scale (£2°), the 1-3% uncertainty in the individual beam transfer 
functions is the largest source of uncertainty, while for multipole moments greater than ~600, random white noise 
from the instrument is the largest source of uncertainty. 


3. THE DATA 

Figure 2 shows the cross-power spectra obtained from all 28 combinations of the 8 differencing assemblies Q1 
through W4 using the quadratic estimator described in Appendix A.l. These spectra have been evaluated with the Kp2 
sky cut described in Bennett et al. (2003c). The spectra are color coded by effective frequency, \J v l v l \ where v [ is the 
frequency of differencing assembly i. The low frequency (4 1 GHz) data are shown in red, the high frequency (94 GHz) 
data in blue, with intermediate frequencies following the colors of the rainbow. The top panel shows /(/+ l)C//27r in 
/xK 2 , while the bottom panel plots the ratio of each channel to our final combined spectrum presented in §5. The top 
panel shows a very robust measurement of the first acoustic peak with a maximum near / ~ 220 and a shape that is 
consistent with the predictions of adiabatic fluctuation models. There is also a clear indication of the rise to a second 
peak at l ~ 540. See Page et al. (2003b) for an analysis and discussion of the peaks and troughs in the first-year WMAP 
power spectrum. 

The red data in the top panel show very clearly that the low frequency data are contaminated by diffuse Galactic 
emission at low l and by point sources at higher l. The higher frequency data show less apparent contamination, 
consistent with the foreground emission being dominated by radio emission, rather than thermal dust emission, as 
expected in this frequency range. 


3.1. Galactic and Extragalactic Foregrounds 


Bennett et al. (2003c) present a detailed model of the Galactic foreground emission based on a Maximum Entropy 
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analysis of all 5 WMAP frequency bands, in combination with external tracer templates. They demonstrate that the 
emission is well modeled by three distinct emission components. 1) Synchrotron emission from cosmic ray electrons, 
with a steeply falling spectrum in the WMAP frequency range: Ta(v) a v& with /? < -3, steepening with increasing 
frequency. 2) Free-free emission from the ionized interstellar medium that is well traced by Ha emission in regions 
where the dust extinction is low. 3) Thermal emission from interstellar dust grains with an emissivity index j} ~ 2.2. 
The model has a Galactic signal minimum between V and W band. 

In principal we could subtract the above model from each WMAP channel and recompute the power spectrum. 
However, since the model is based on WMAP data that have been smoothed to an angular resolution of 1 .°0, the 
resulting maps would have complicated noise properties. For the purposes of power spectrum analysis, we adopt a 
more straightforward approach based on fitting foreground tracer templates to the Q, V, and W band data. The details 
of this procedure, the resulting fit coefficients, and a comparison of the fits to the Maximum Entropy model are given 
in Bennett et al. (2003c). They estimate that the template model removes ~85% of the foreground emission in Q, V, 
and W bands and that the remaining emission constitutes less than ~2% of the CMB variance (up to l = 200) in Q 
band, and < 1 % of the CMB variance in V and W bands. 

The contribution from extragalactic radio sources has been analyzed in three separate ways. Bennett et al. (2003c) 
directly fit for sources in the sky maps. The result of this analysis is that we have identified 208 sources in the WMAP 
data with sufficient signal to noise ratio to pass the detection criterion (we estimate that ~5 of these are likely to 
be spurious). The derived source count law is consistent with the following model for the power spectrum of the 
unresolved sources 



with A = 0.015 }jlK 2 sr (measured in thermodynamic temperature), ft = -2.0, and = 45 GHz. Komatsu et al. (2003) 
evaluate the bispectrum of the WMAP data and are able to fit a non-Gaussian source component based on a particular 
configuration of the bispectrum. They find the same source model, equation (9), fits the bispectrum data. For the 
remainder of this section, we adopt this model for correcting the cross-power spectra. At Q band (41 GHz) the 
correction to /(/+ l)C//27r is 868 and 3468 fi K 2 at /=500 and 1000, respectively. At W band (94 GHz), the correction 
is only 31 and 126 /iK 2 at the same l values. For comparison, the CMB power in this / range is ~ 2000 //K 2 . Later, 
when we derive a final combined spectrum from the multi-frequency data, we adopt equation (9) as a model with A as 
a free parameter. We simultaneously fit for a combined CMB spectrum and source amplitude and marginalize over the 
residual uncertainty in A. The best-fit source amplitude from this process is consistent with the other two methods. 

Figure 3 shows the cross-power spectra obtained from the same 28 combinations as in Figure 2, this time with 
the Galactic template model and source model subtracted. The bottom panel of the Figure shows the ratio of the 28 
channels to the combined spectrum obtained in §5. The 28 cross-power spectra are consistent with each other at the 
5 to 20% level over the / range 2 - 500. Similar scatter is seen in Monte Carlo simulations of an ensemble of 28 
cross-power spectra with WMAFs beam and noise properties. The only significant deviation lies in the Q band data at 
low l which is ~ 1 0% below the higher frequency bands at l < 20. This is consistent with the accuracy estimated above 
for the Galactic template model, see also Figure 1 1 of Bennett et al. (2003c) for images of the maps after Galactic 
template subtraction. Since the WMAP data are not noise limited at low /, we use only V and W band data in the final 
combined spectrum for l < 100. 

The subtraction of the source model, equation (9), brings the Q band spectrum into good agreement with the other 
cross-power spectra up to l ~ 500. At higher /, the Q band data contributes very little to the final combined spectrum 
because the (normalized) Q band window function has dropped to less than 5% (Page et al. 2003a). As discussed in 
§5, we marginalize over the source amplitude uncertainty, SA , when obtaining the final power spectrum estimate and 
associated covariance matrix. Thus the uncertainty is also accounted for in subsequent cosmological parameter fits 
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( Verde et al. 2003; Spergel et al. 2003; Peiris et al. 2003). 

Figure 4 shows a close-up of the 28 cross-power spectra in Figure 3 up to l = 100. The top panel shows the raw 
(unbinned) data which has correlations of <2% between neighboring points is this l range (see §4). These data are 
strikingly consistent with each other and support the conclusion that systematic errors at low l are insignificant. To 
assess the level of scatter that does exist between the spectra, we have generated a Monte Carlo simulation in which we 
compute the rms scatter among the 28 spectra at each L relative to the measured power. The bottom panel of Figure 4 
shows the results of this simulation, averaged over 1000 realizations, compared to the relative rms scatter in the data. 
The agreement is excellent, indicating that the uncertainty in the measured power spectrum in this l range is a few 
percent and is consistent with a combination of instrument noise and mode coupling due to the 15% sky cut. 

Another striking feature is the low amplitude of the observed quadrupole, and the sharp rise in power, almost 
linear in /, to / = 5. Bennett et al. (2003b) quote a value for the rms quadrupole amplitude, Q rms = -y/(5/47r)C2 = 8 ± 2 
H K, where the uncertainty is largely due to Galactic model uncertainty. This is consistent with the amplitude measured 
by the COBE-DMR experiment, Qrms = 10^ /xK (Bennett et al. 1996). The fast, nearly linear rise to / = 5 produces 
an angular correlation function with essentially no power on angular scales > 60°, again in excellent agreement with 
the COBE-DMR correlation function (Bennett et al. 2003b; Hinshaw et al. 1996). In the context of a standard ACDM 
model, the probability of observing this little power on scales greater than 60° is ~ 2 x 10 -3 (Bennett et al. 2003b). 


4. THE FULL COVARIANCE MATRIX 

In §5 we derive our best estimate of the angular power spectrum by optimally combining the 28 cross-power 
spectra discussed above. The procedure for combining spectra requires the full covariance matrix of the individual 
cross-power spectra - in this section we outline the salient features of this matrix. There are six principal sources of 
variance for the measured spectra, Cj: cosmic variance, instrument noise, mode coupling due to the foreground mask, 
point source subtraction errors, uncertainty in the beam window functions, and an overall calibration uncertainty. 
We ignore uncertainties in the diffuse foreground correction since they are everywhere sub-dominant to the cosmic 
variance uncertainty (see §3.1). 

We may write the covariance matrix as 

(Efu4< = ([q-(C,+ASVj][Cj, ~(C r +ASV;, ]) , (10) 

where the angle brackets represent an ensemble average, C/ is the true underlying power spectrum, w) is the window 
function of spectrum i, and AS' is the point source contribution to spectrum i. Here we have defined a point source 
spectral function, S', as 



where v x , v$, and ft are as defined after equation (9). Note that (SfuiOjj, is symmetric in both (//') and (ij). 

In the process of forming the combined spectrum we will estimate a best-fit point source amplitude, A, and 
subtract the corresponding source contribution from each spectrum i. We thus rewrite Efuii as 

(Sfuii)li = {[Cj - (G+SASV; ][Cj, - (C// + 5AS>)J r ]) , (12) 

where Cj = Cj-ASV, is the source-subtracted spectrum, and 6A = A- A is the residual source amplitude, which we 
will marginalize over as a nuisance parameter. 
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We may expand the covariance matrix as 


— ^cv + ^mask + ^src + 


(13) 


We discuss each of these contributions in more detail below. 

Cosmic Variance , Noise , and Mode Coupling. — The first two terms, E cv + E mas k, incorporate the combined 
uncertainty due to cosmic variance, instrument noise, and mode coupling due to the foreground mask, 

(Scv + Smask)?;, = ([C?-Qw}][Cf, ~C V *{,]), _ (14) 

where w\ is fixed at its measured value. We have split this contribution into two pieces to mimic the procedure we 
actually use to compute the covariance matrix. As outlined in §5.1, we start with E cv , then incorporate the effects of 
point source error and window function error. We do not add the effects of mode coupling, Emask* unt il the very end 
of the computation. We consider this term in more detail in Appendix D. 

Point Source Subtraction Errors . — The third term, E src , is due to uncertainty in the point source amplitude 
determination 

(S src ))j, = 5V;<7 s 2 rc 45j, (15) 

where cr 2 c = ((A -A) 2 ) is the variance in the best-fit amplitude A, and we assume that the frequency dependence, S\ 
is perfectly known. In practice, we do not explicitly evaluate E src as given above, rather we employ a method based on 
marginalizing a Gaussian likelihood function, £(CJ, A), over a nuisance parameter A. This process, which is discussed 
in Appendix B, effectively yields E cv + E src . 

Window Function and Calibration Errors . — The fourth term, E W i n , is due to uncertainty in the beam window 
function, vv}. This term arises from fluctuations in the window function which cause the measured spectrum, Cj to differ 
from our estimate of the convolved spectrum, C/wJ, where w) is the estimated window function. This contribution has 
the form 

(E win )« = C,(Aw‘.a4)c ; -. (16) 

Recall from §2 that w) = b\b{ pj where b\ is the beam transfer function for DA i and p } is the pixel window function. 
Define u\ = A b\jb\ to be the fractional error in b\ , then to first order in u\ we have 

^AwJ * Aw},^ = w) ^ul+u^uj, +uj uj, + «/ uj, ^ wj, . (17) 

For WMAPthe uncertainty in b\ is uncorrelated between DAs, thus the above expression reduces to 

^Aw} * Aw},^ =wj ^uji’iSij+Sij^+Ti^iitiSjij+Siij^ wj, , (IB) 

where we have defined 

and where w is the covariance matrix of the beam transfer function for DA i given by Page et al. (2003a). When 
generating the combined spectrum in the next section, we add E W i n to the above contributions, giving E cv + E src + E W in- 

The WMAP absolute calibration uncertainty is 0.5%. We do not explicitly incorporate this contribution in the 
covariance matrix. Instead, we propagate a 0.5% uncertainty into the normalization of the final power spectrum 
amplitude (Spergel et al. 2003). 
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5. THE COMBINED POWER SPECTRUM 

In §3 we demonstrated that the three high frequency bands of WMAP data produced consistent estimates of 
the angular power spectrum, after a modest correction for diffuse Galactic emission and extragalactic point sources. 
It is therefore justifiable to combine these data into a single “optimal” estimate of the angular power spectrum of 
the CMB. In this section, we provide an overview of two methods we use to generate a single combined spectrum. 
The first is a multi-step process that simultaneously fits the 28 cross-power spectra presented above to a single CMB 
power spectrum and a point source model, equation (9), while correctly propagating beam and residual point source 
uncertainties through to a final Fisher matrix. This spectrum constitutes our best estimate of the CMB power spectrum 
from the first-year WMAP data. The second spectrum, which serves as a cross check of the first, is based on forming a 
single co-added sky map from the Q1 through W4 maps, and using the quadratic estimator with noise bias subtraction. 
We compare the two spectra in §5 .3. 


5.1. Method I - Optimal Combination of Cross-Power Spectra 

Since this method is relatively complicated, we outline the basic procedure here and relegate the details to Ap- 
pendices, as indicated. We present the result in §5.3. The steps are as follows. 

1. Subtract best-fit Galactic foreground templates from each of the maps Q1 through W4, using the coefficients 
given in Table 3 of Bennett et al. (2003c). 

2. Evaluate the 28 cross-power spectra from the maps Q1 through W4, where each spectrum has been evaluated 
using the quadratic estimator of Appendix A.l with the weighting scheme defined in Appendix A. 1.2. 

3. Collect the noise bias estimate, (Nj), for each DA from §2.2. These estimates are used in the calculation of the 
covariance matrix for the combined spectrum, and in setting the relative weight of each cross-power spectrum 
in the final combined spectrum. 

4. Apply the procedure presented in Appendix B.l to obtain an estimate for the point source amplitude. The value 
obtained is A = 0.0155 ± 0.0017, roughly independent of / max in equation (B6). This value is in good agreement 
with an estimate based on the bispectrum (Komatsu et al. 2003), and on an extrapolation of point source counts 
(Bennett et al. 2003b). Subtract the point source contribution from the cross-power spectra: C] = Cj -A^wj. 

5. Compute an approximate form of the full covariance matrix discussed in §4, E^u. The procedure we use pro- 
duces a covariance matrix which includes cosmic variance, instrument noise, source subtraction uncertainties, 
and window function uncertainties. At this stage of the process, it does not yet include the effects of mode 
coupling. More details are given in §4 and Appendix B. 

6. Invert the approximate covariance matrix for use in computing the optimal spectrum. This is the most computa- 
tionally intensive step in the process. 

7. Compute the final combined spectrum from the 28 CJ as per the procedure given in Appendix C. In particular, 
assume a fiducial cosmological model (as specified in the Appendix), and use equation (C6), with E^ w = 
Ejjijj lv . This produces a final spectrum which is very nearly optimal. 

Note: for l < 100 we use a surrogate procedure for computing the combined spectrum. In order to minimize 
Galactic foreground contamination, we use only V and W band data. Moreover, because the statistics of the Cy 
are mildly non-Gaussian, and because point source contamination and window function uncertainties are small, 
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the “optimal” machinery developed above is unnecessary. Rather we simply form a weighted average spectrum 
from the V and W band Cj. 

8. Compute the approximate inverse-covariance matrix, Q , for the combined spectrum using equation (C7). This 
matrix is approximate in two ways. 1) It does not yet incorporate the effects of mode coupling - this is added 
below. 2) It has been evaluated for a fiducial cosmological model, Q(Cf d ), while in a likelihood application, we 
need to evaluate Q{Cf) for an arbitrary model, Cf - we add this next. 

9. Introduce the dependence on cosmological model into Q as follows. Invert Q to obtain the approximate co- 
variance matrix of the combined spectrum, E. The off-diagonal terms of E are small and weakly dependent on 
cosmological model, so we expand E as 

E ut =Didii'+€u /, (20) 

where €/// encodes the mode coupling due to window function and source subtraction uncertainties. This relation 
defines €///, which we take to be zero on the diagonal. Dj is dominated by cosmic variance and noise; in order 
to separate the two contributions, we introduce an ansatz of the form 

<2i> 

where Cf d is the fiducial model used to generate the combined spectrum (see Appendix C), and Nf is the 
effective noise in the combined spectrum, which is defined by this equation. The dependence on cosmological 
model is introduced in the covariance matrix by re-computing £>/ with Cf d Cf\ leaving Nf n fixed. 

10. Estimate the coupling induced by the foreground mask as described in Appendix D. The effect of the mask on 
the off-diagonal elements of the Fisher (inverse-covariance) matrix can be written as 

(E -^) w =r w (D' l D' v T x/ \ (22) 

where D\ = This expression parameterizes, and defines, the off-diagonal mode coupling in the form of a 

correlation matrix, r w . Note that F w , defined in Appendix D.2, is related to r iV by Fu< = + ru* . 

11. The final Fisher (curvature) matrix, Qu^ is obtained using 

Qw =(D' I r 1 5„.-e w (DW') l / 2 + r w (DWr l/ \ (23) 

We further calibrate D\ with Monte Carlo simulations. This calibration process, and a description of how the 
curvature matrix is used in a maximum likelihood determination of cosmological parameters, is given in §2 of Verde 
et al. (2003). As part of the first-year data release, we provide a Fortran 90 subroutine to evaluate the likelihood of a 
given cosmological model, C{ h , given the WMAPdata (supplied in the routine). The code also optionally returns the 
Fisher (inverse-covariance) matrix for the combined spectrum. 


5.2. Method II - Combined Sky Map 


Our second approach is to form a single co-added map from the Q 1 through W4 maps. 


T = 




(24) 
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where T{ is the sky map for DA i with the best-fit Galactic template model subtracted, and a^ is the noise per 
observation for DA i, given by Bennett et al. (2003b). We evaluate the power spectrum of this map on the Kp2 cut sky 
using the quadratic estimator in Appendix A.l. An effective noise model is obtained by using the same approach as 
described in §2.2: we generate co-added noise maps from the library of end-to-end simulations, evaluate their average 
spectra, then fit a noise model. The noise bias model is then subtracted from the power spectrum of the combined 
temperature map. We have performed this analysis with 3 distinct pixel weighting schemes (see Appendix A. 1.2) and 
3 corresponding noise models. The results are shown in Figure 5 where it is seen that the three cases are virtually 
indistinguishable. 

In effect, this analysis uses both the auto- and cross-power spectra. We view this as a useful check of the more 
sophisticated procedure described in §5.1, but we do not rely on it for a final result. Uncertainties in the noise model 
only effect the fourth moment of the cross-power spectra, but they effect the second moment of the auto-power spectra 
and potentially bias the final result. The ~6% sensitivity advantage gained by including auto-power spectra was not 
deemed worth the effort required to guarantee that the final result was not biased. 


5.3. Comparison of Results 

Figure 6 compares the power spectra obtained from methods I and II above. The combined cross-power spectrum 
from §5.1 is shown in black, the auto-power spectrum obtained in §5.2 from the co-added map is shown in grey. The 
two methods agree extremely well with the only notable deviation being at the highest l range probed by the first-year 
data. This is the regime where the auto-power spectrum will be most sensitive to the noise bias subtraction. As can be 
seen in the error estimates shown in Figure 8, the deviation is less than la. 

A separate test of robustness is to compute the angular power spectrum in separate regions of the sky to see if 
the spectrum changes. We have computed the power spectrum in two subsets of the sky - the ecliptic poles, and the 
ecliptic plane, using the quadratic estimator with the combined sky map. The results are shown in Figure 7 where the 
pole data is shown in grey and the plane data in black. The two spectra are very consistent overall, but some of the 
features that appear in the combined spectrum, such as the “peak” at / ~ 40 and the “bite” at / ~ 210, are not robust to 
this test, thus we consider these features to be of marginal significance. There is also no evidence that beam eliipticity, 
which would be more manifest in the plane than in the poles, systematically biases the spectrum. This is consistent 
with estimates of the effect given by Page et al. (2003a). 


6. DISCUSSION 

Our best estimate of the angular power spectrum of the CMB is shown in Figure 8. Also shown is the best-fit 
ACDM model from Spergel et al. (2003) which is based on a fit to the this spectrum plus a compilation of additional 
CMB and large-scale structure data. The WMAPdata points are plotted with measurement errors based on the diagonal 
elements of the Fisher matrix presented in Appendix D. The cosmic variance errors, which include the effects of the 
sky cut, are plotted as a la band around the best-fit model. As discussed in Spergel et al. (2003), the model is an 
excellent fit to the data. The combined spectrum provides a definitive measurement of the CMB power spectrum 
with uncertainties limited by cosmic variance up to / ^ 350. The spectrum clearly exhibits a first acoustic peak at 
l = 220.1 ± 0.8 and a second acoustic peak at l = 546 ± 10. Page et al. (2003b) present an analysis and interpretation 
of the peaks and troughs in the first-year WMAP power spectrum. 

Figure 9 compares the first-year WMAP spectrum to a compilation of recent balloon and ground-based measure- 
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ments. In order to make this Figure meaningful, we plot the best-fit model spectrum to represent the WMAP results. 
The data points are plotted with errors that include both measurement uncertainty and cosmic variance, so no error 
band is included with the model curve. (Since individual groups report band power measurement with different band- 
widths, it is not possible to represent a single cosmic variance band that applies to all data sets.) The model spectrum 
fit to WMAP agrees very well with the ensemble of previous observations. 

Wang et al. (2002a) have recently distilled a CMB power spectrum from an optimal combination of the extent 
pre- WMAP data. In Figure 10 we plot their derived band power points alongside the WMAP data. To make this 
comparison meaningful, we plot the WMAP data with cosmic variance plus measurement errors and omit the error 
band from the model spectrum. The distilled spectrum is notably lower than the WMAP data in the vicinity of the 
first acoustic peak. In a previous version of this work (Wang et al. 2002b) the authors noted that the first peak of their 
combined spectrum was lower than a significant fraction of their input data. They attribute this to their formalism 
allowing for a renormalization of individual experiments within their respective calibration uncertainties. Figure 1 in 
Bennett et al. (2003a) presents a similarly distilled spectrum from the data extent in late 2001 and found a first peak 
amplitude that was more intuitively consistent with the bulk of the input data, and which is now seen to be consistent 
with the WMAP power spectrum. 

Figure 1 1 shows the WMAP combined power spectrum compared to the locus of predicted spectra, in red, based 
on a joint analysis of pre- WMAP CMB data and 2dFGRS large-scale structure data (Percival et al. 2002). As in 
Figure 8, the WMAP data are plotted with measurement uncertainties, and the best-fit ACDM model (Spergel et al. 
2003) is plotted with ala cosmic variance error band. Percival et al. (2002) predict the location of the first peak should 
occur at l = 221.8 ±2.4, which is quite consistent with the value reported by Page et al. (2003b) of / = 220.1 ±0.8. The 
height of the first peak was predicted to be in the range 4920 ± 170 /zK 2 , while Page et al. (2003b) report a measured 
height of 5580 ± 75 fjK 2 , about 13% higher. Unlike the position, the amplitude of the first peak has a complicated 
dependence on cosmological parameters. Percival et al. (2002) report best-fit parameters for a ACDM model that are 
mostly consistent with those reported by Spergel et al. (2003) for the same class of models. The only mildly disparate 
comparison lies in the combination of normalization, as, and optical depth, r. Percival et al. (2002) report the product 
v%e~ T = 0.72 ± 0.03 ± 0.02, where the first error is a “theory” error and the second is measurement error. While 
Spergel et al. (2003) does not report a maximum likelihood range for this explicit parameter combination, the product 
of their maximum likelihood values for ag and r yields a%e~ T = 0.74, which is consistent with Percival et al. (2002), 
but would make the first peak a few percent higher. Small differences in £l m h 2 y and n s , may also contribute to 

the difference. 


7. CONCLUSIONS 

We present measurements of the angular power spectrum of the cosmic microwave background from the first- 
year WMAP data. The eight high-frequency sky maps from DAs Q1 through W4 were used to estimate 28 cross- 
power spectra, which are largely independent of the noise properties of the experiment. These data were tested for 
consistency in §3, then used in §5 as input to a final combined spectrum, discussed in §6. The procedure for estimating 
the uncertainties in the final combined spectrum were discussed in §4 and in numerous Appendices. 

The combined spectrum provides a definitive measurement of the CMB power spectrum, with uncertainties lim- 
ited by cosmic variance up to / ~ 350, and a signal to noise per mode > 1 up to l ~ 650. The spectrum clearly exhibits 
a first acoustic peak at l = 220.1 ±0.8 and a second acoustic peak at / = 546 ± 10. Page et al. (2003b) present an 
analysis and interpretation of the peaks and troughs in the first-year WMAP power spectrum. Spergel et al. (2003), 
Verde et al (2003), and Peiris et al. (2003) analyze the combined spectrum in the context of cosmological models. 
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They conclude that the data provide strong support adiabatic initial conditions, and they give precise measurements of 
a number of cosmological parameters. Kogut et al. (2003) analyze the correlation between WMAPs temperature and 
polarization signals, the Cf E spectrum, and present evidence for a relatively high optical depth, and an early period 
of cosmic reionization. Among other things, this result implies that the temperature power spectrum is suppressed by 
~30% on degree angular scales, due to secondary scattering. 

A variety of first-year WMAP data products are being made available by NASA’s new Legacy Archive for Mi- 
crowave Background Data Analysis (LAMBDA). In addition to the sky maps and calibrated time-ordered data, we 
are providing the 28 cross power spectra used in this paper (with diffuse foregrounds subtracted), the combined spec- 
trum from §5.1, and a Fortran 90 subroutine to compute the likelihood of a given cosmological model, (the code will 
also optionally return the Fisher (inverse-covariance) matrix for the combined spectrum.) The LAMBDA URL is 
http: //lambda . gsfc . nasa . gov/. 
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of this manuscript. LV is supported by NASA through Chandra Fellowship PF2-30022 issued by the Chandra X-ray 
Observatory center, which is operated by the Smithsonian Astrophysical Observatory for an on behalf of NASA under 
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A, POWER SPECTRUM ESTIMATION METHODS 

For the analysis of WMAF s first-year data, we have chosen two distinct methods for inferring the power spec- 
trum. The first is a fast and accurate quadratic method for estimating the power spectrum of a partial sky map (Hivon 
et al. 2002). We summarize the basic approach here, highlighting the aspects of the method that are especially perti- 
nent to WMAP, and refer the reader to Hivon et al. (2002) for details. The second is a maximum likelihood method 
that provides an independent estimate of the spectrum measured by WMAP (Oh et al. 1999). 


This preprint was prepared with the AAS I^TgX macros v5.0. 
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We have applied both of these methods to the WMAP data. The results are shown in Figure 12, which shows 
spectra estimated from the V band map, up to l = 200, for the two methods. The maximum likelihood estimate has 
slightly lower uncertainties at low /, because the method optimally weights the data with a pixel-pixel covariance 
C = S+N « S where S is the covariance of the CMB signal and N is the covariance of the noise (see Appendix A.2). 
Our quadratic estimator uses uniform pixel weights at low l (see Appendix A. 1.2) though it is clear from the Figure 
that the difference is not significant. At high /, where the WMAP data are noise dominated, the two estimators give 
essentially identical results because they effectively weight the data in the same way. 

To obtain our “best” estimate of the WMAP power spectrum, we adopt the quadratic estimator because it can be 
readily applied to pairs of WMAP radiometers in a way that is nearly independent of the properties of the instrument 
noise. In §5 we discuss our methodology for combining spectra from pairs of radiometers and present the final 
combined spectrum. 


A.l. Quadratic Estimation 

Hivon et al. (2002) start with the full-sky estimator in equation (2), add a position dependent weight, W(n), and 
set W to zero in the regions where the sky is contaminated. In other regions, W is chosen to optimize the sensitivity of 
the estimator. A temperature map AT (n) on which a weight W(n) is applied can be decomposed in spherical harmonic 
coefficients as 

aim = J dn n AT(n)W(n)Y,*J n ) (Al) 

« %^2at(p)W(p)y;jp), (A2> 

P 

where the integral over the sky is approximated by a discrete sum over map pixels, each of which subtends solid angle 
Hivon et al. (2002) then define the “pseudo power spectrum” Q as 

i / 

c ' = 27Ti£M 2 - (A3) 


The pseudo power spectrum C/, given by the weighted spherical harmonic transform of a map, is clearly different 
from the full sky angular spectrum, C ; sky , but the ensemble averages of the two spectra can be related by 


<q)=£G„,<cf) (A4) 

V 

where Gu> describes the mode coupling resulting from the weight function W(n) (Hauser & Peebles 1973). Hivon 
et al. (2002) give the following expression for the coupling matrix, which depends only on the geometry of the weight 
function W(n) 

o o/' (A5) 

where the final term in parentheses is the Wigner 3- j symbol, and W/ is the angular power spectrum of the weight 
function 

W = 5T T £ k 'I ! 


(A6) 
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where 

w lm = j <m a W(n)Y,* m (n). (A7) 

Upon inverting the coupling matrix G/// and making the identification (Cf*) = C/, we obtain the following esti- 
mator of the power spectrum 

C?* = Y, G 7i'Ci'- (A8) 

v 

The computation of equation (A2) for each (/,m) up to l = / max would scale as Npix^ax if performed on an arbitrary 
pixelization of the sphere, where N p j X is the number of sky map pixels. However, for a pixelization scheme with 
iso-latitude pixel centers, fast FFT methods may be employed to speed up the evaluation, so it scales like 
(Muciaccia et al. 1997). The WMAP sky maps have been produced using the HEALPix layout (Gorski et al. 1998) 
which supports such fast spherical harmonic transforms. In particular, the HEALPix routine map2alm evaluates 
equation (A2). 


ALL Auto- and Cross-Power Spectra from the WMAP Data 

For a multi-channel experiment like WMAP it is quite powerful to evaluate the power spectra from different maps 
and compare results. In particular, the quadratic estimator described above may be used on 1 or 2 maps at a time by 
generalizing the expression for the pseudo power spectrum equation (A3) as 

(A9> 

m=-l 

where a\ m refers to the transform of map i and d [ r * refers to map j , which needn’t be the same as map i. As discussed 
in §2, if i f- j and the noise in the two maps is uncorrelated, the estimator equation (A9) provides an unbiased estimate 
of the underlying power spectrum. 

We have tested the auto- and cross-power estimator extensively with Monte Carlo simulations of the first-year 
WMAP data. Selected results from this testing are shown in Figure 13. The auto- and cross-power spectra that obtain 
from the flight data are presented in detail in §3. The cross-power spectra form the basis for our final combined 
spectrum, presented in §5. 


A. 7. 2. Choice of Weighting 

We seek a weighting scheme that mimics the maximum likelihood estimation (Appendix A.2), which effectively 
weights the data by the full inverse covariance matrix C -1 = (S + N) _1 . For the combined spectrum presented in §5, we 
use three distinct weighting functions in three separate / ranges. 

1 . For / < 200 we give equal weight to all un-cut pixels, 

W(p) = M(p) (A10) 

where M(p) is the Kp2 sky mask, defined by Bennett et al. (2003c). It takes values of 0 within the mask and 1 
otherwise. 


2. For / > 450 we use inverse-noise weighting, 


W(p) = M(p)N obs (p) 


where N obs (p) is the number of observations of pixel p. 

3. For 200 < / < 450 we use a transitional weighting, 

W(p). ^ 

P l/(N obs ) + l/N obs (p) 

where {N 0 b S ) is the mean of W 0 bs evaluated over the cut sky. 


Verde et al. (2003) discuss the choice of weighting in more detail. 


(All) 


(A 12) 


A.2. Maximum Likelihood Estimation 


This Appendix provides a summary of the maximum likeli hood approach to power spectrum estimation originally 
presented by Oh et al. (1999). If the temperature fluctuations are Gaussian, and the a priori probability of a given set 
of cosmological parameters is uniform, then the power spectrum may be estimated by maximizing the multi-variate 
Gaussian likelihood function 


£(C;|m) 


exp(-jm r C ’m) 
%/detC 


(A13) 


where m is a data vector (see below) and C is the covariance matrix of the data, which has contributions from both 
the signal and the instrument noise, C = S+N. We can work in whatever basis is most convenient. In the pixel basis 
the data are the sky map pixel temperatures, and in the spherical harmonic basis the data are the ai m coefficients of the 
map. In the former basis the noise covariance is nearly diagonal, while in the latter, the signal covariance is. 


Pixel basis: Spherical harmonic basis: 

m — y Ti m —> oi m 

s -4 £/ ^C,fi(cos0 y ) S -»• diag(C 2 ,C 2 ,...C 3 ,C 3 ,...) 

N -*• afSij N -+ N ( i m)i i m y (see below). 

For WMAP, the length of the data vector, A/data, is 2,672,36 1 , the number of 7' sky map pixels (HEALPix V S ide = 512) 
that survive the Galaxy cut, so it is necessary to find methods for evaluating C that do not require a full inversion of 
the covariance matrix C, which requires O (A^) operations. We use an iterative method for evaluating the likelihood 
that exploits the ability to find an approximate inverse C" 1 . The most important features in the data that make this 
possible are that WMAP observes the full sky and the Galaxy cut is predominantly azimuthally symmetric in Galactic 
coordinates (Bennett et al. 2003b). Of secondary importance for this pre-conditioner is that WMAF s noise per pixel 
does not vary strongly across the sky (Bennett et al. 2003a). We discuss the pre-conditioner in more detail below. 

Defining / = ~2ln£ and P/ = J^, we maximize the likelihood by solving 

^ = 0 = m r C' I P ; C- , m+tr(C- , P,) (A14) 

oQ 

using a Newton-Raphson root finding method that generates an iterative estimate of the angular power spectrum 
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at each step. Here Fu> is the Fisher matrix 

F w =-(S 2 ln C/dQdCr) = ^trCCT 1 P/CT 1 P/0- (A16) 

To implement the solution in equation (A 15) we need a fast way to evaluate the following three components of 
'Li' F w§c l 7 : 

1. m^C 'P/C 'm 

2. tr(C-'P,) 

3. tr(C -1 P/C -1 P//). 


We use the spherical harmonic basis in which the data vector consists of the ai m coefficients of the map obtained by 
least squares fitting on the cut sky. The signal covariance is diagonal in this basis, while the noise matrix is obtained 
from the normal equations for the ai m fit 

^ a (lmY = ?(/«) (A17) 

(imy 

where 


N(im)(i m y 

yam) 


E 

i 

E 


T(/m) (ImY (Ai) 

V] 

rrT’ * 


(A 18) 
(A19) 


The sums are over all sky map pixels that survive the Galaxy cut, and we have assumed that the noise is uncorrelated 
from pixel to pixel. 


A.2.7. Evaluation of C 1 m 


The term C 1 m appears repeatedly in the evaluation of equation (A 15). We compute this by solving Cz = 
(S + N) z = m for z. A more numerically tractable system is obtained by multiplying both sides by S 2 N" 1 , so 


(I+S^ N -1 S^)S 2 2 = S 2 N _1 m = y 


(A20) 


where y is the spherical harmonic transform of the map, defined in equation (A 19). Note that y can be quickly 
computed in any pixelization scheme that has iso-latitude pixel centers with fixed longitude spacing. We then solve 
equation (A20) using an iterative conjugate gradient method with a pre-conditioner for the matrix A = (I+S 3 N" 1 ). 

We find the following block-diagonal form of A to be a good starting point 


A = 


I+S^N-'S 2 

0 


diag(I+S^N-'S*) 


(A21) 


where N" 1 is a block-diagonal approximation of the noise matrix discussed below. The lower-right block of A occupies 
the high / portion of the matrix where the signal to noise ratio S- N~ ! is low, so a diagonal approximation is 
adequate. The upper-left block occupies the low / portion of the matrix where the signal dominates the noise, so we 
need a better estimate of N” 1 . In practice we find this split works well at / = 512 for the WMAP noise levels. As for 
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the approximate form of N“ l , defined in equation (A19), note that the dominant off-diagonal contributions arise from 
the azimuthally symmetric Galaxy cut, which couples different / modes, but not m modes. Thus N -1 is approximately 
block diagonal, with perturbations induced by the non-uniform sky coverage of WMAP . We therefore use a block 
diagonal form of N~ } as the pre-conditioner, 

= $rnm' - (A22) 

Using the pre-conditioner equation (A21) we find that the conjugate gradient solution of equation (A20) converges in 
approximately six iterations and requires only cpu-minutes of processing on an SGI Origin 2000. 


B. POINT SOURCE SUBTRACTION 


In this Appendix we describe the procedure we use to estimate and subtract the point source contribution directly 
from the multi-frequency cross-power spectra. We then show how we incorporate the source model uncertainty into 
the covariance matrix of the source-subtracted spectra by marginalizing a Gaussian likelihood function over the source 
model amplitude parameter. 


This marks the starting point of the multi-frequency analysis which will lead to the combined power spectrum, 
discussed in §5.1. In order to generate the combined spectrum, we need the full covariance matrix of the cross-power 
spectra (see §4). Our approach to generating the full covariance is to start with the ideal, full-sky form, which only 
includes cosmic variance and instrument noise, then we incorporate additional effects step by step, as outlined in §5.1 
and in these Appendices. For an ideal experiment with full sky coverage, no point source contamination, and no beam 
uncertainty the covariance matrix is 


v" JJ - 
11 ( 2 / + 1 ) 


h w 1 / + n l nj Si j'j ; +n l n J 5? + (c^w 1 / +n l n J 8\y j 7 + (B 1) 


where denotes the Kronecker symbol, w\ } = b\b\ pf is the window function, and n l n l 5# = NT 


B.l. Estimating the Point Source Amplitude 

We start by assuming a Gaussian likelihood for the sky model, given the WMAP data 

-21n£(A,C ; |C*) = £ [Cj-(C,+ASy}] (E- ] )j j [Cj-(C,+ASVj], (B2) 

y/ 

where Cj is cross-power spectrum i, w ) is the window function for spectrum i, C/ is the true CMB power spectrum, and 
AS* is the source model defined in equations (9) and (11). Here we assume the diagonal form of the covariance matrix 
in equation (Bl). 

To determine the best-fit source amplitude. A, we marginalize this likelihood over the CMB spectrum, C/. First 
we expand equation (B2) as 

-2 In C = 

l u 

- Y, 2C ; Y, (Cj -ASX XE- 1 )jx 
/ ij 

+ Y C ‘Y 

l u 


(B3) 
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which is a quadratic form in C/, C a f[/ exp [~j(aCf - 2 bQ+c)], with 

>j 

b = ^(q-AS'w’XE- 1 )^, 

U 

c = ^(C7-A5V;XS*‘); j (Cf-A5 j wi). (B4) 

U 


We wish to evaluate the marginalized likelihood, Cc t {A) = f dQC(A,Ci). This is readily evaluated using the substitu- 
tion Q -► Q-b/a , giving Cq oc ri/ ex P[ - i( c ~^ 2 / fl )]» where we drop a multiplicative term proportional to a , which 
is independent of A, The marginalized likelihood function is thus 


-21n£ C/ 


£ £(Cj -ASVjXE - 1 ) ; S (C? -ASVj) 
/ ij 



1 2 

^(Cj-ASVjXE-'W • 

u 


(B5) 


Setting d£c,(A)/dA = 0 gives the most likely point sources amplitude 


where 


H‘ = S l - 


Epj, »f (E- | )|' , wi' 


The standard error on the best fit value of A is 


(B6) 


(B7) 


cr 


2 _ 
src — 


«/ 


-1 


(B8) 


For WMAP , the off-diagonal terms in the covariance matrix are small. Here, neglecting them changes the inferred 
point source amplitude by less than 


B.2. Marginalizing over Point Source Amplitude 

The source subtraction procedure discussed above is uncertain. In this section we incorporate this uncertainty 
into the full covariance matrix for the cross-power spectra. We again assume a Gaussian likelihood function of the 
form 

-21n£ = ^ [C;-(Cf + *ASj)w}](£- , )|{, [4-(C?+M4H] (B9) 

ij W 

where Cj is the source-subtracted cross-power spectrum for DA pair i, obtained above, w\ is the window function for 
spectrum i, CJ h is now a fixed CMB model spectrum, and SA is the residual source amplitude. 
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We can marginalize the likelihood function over the residual point source amplitude as follows. Expand equa- 
tion (B9) as 

— 21n C = Y -qv;, ) 

ij IV 

- 2(6A)Y 

ij IV 

+ (SA) 2 Y (bio) 

ij//' 

which is a quadratic form in <5A, C oc exp[-^(a(5A) 2 -2b(5A)+c)], with 

a = Y s ' w ' 9 

ij iv 

b = ^(q-cfwjxE 

ij //' 

c = £(3-<^x2r 1 )^ (bid 

ij //' 

We wish to evaluate the marginalized likelihood, C A - f Cd(SA). This is readily evaluated using the substitution 
SA -> 5A-b/a, giving C A cx exp \-\{c-b 2 fa)], where we drop a multiplicative term proportional to a, which is 
independent of q. The marginalized likelihood function is thus 

-21n£ A = ^<q-c; h vw;)(s- , )5j,(4-q?wj,) (B12) 

ij IV 

- -^(q-q h w;xs-')|j,s j 4 Y <bi3> 

a ij iv vy vv " 


This expression may be recast in the form 

-21n£ A = Y (q-q h w|) (F &T % (4 -C?4 ), (B14) 

ij IV 


where (F src )jj, is 

(F src );j,=(E -' )%- ±B%, (B15) 

with 

B% = Y (S -1 )ii'»/wf,,(E- 1 )]j',,,,5 j '4,,. (B16) 

j/j/ [ 1 * 1 *" 

The superscript “sre” indicates that the Fisher matrix so obtained includes point source subtraction uncertainty, in 
addition to whatever effects have been included in (E -1 )/j, already, in this case only cosmic variance and noise. Note 
that equation (B15) neglects a term proportional to lndeta which has a weak dependence on cosmological parame- 
ters. In the actual calculation, as in the previous section, we assume the diagonal form of the covariance matrix in 
equation (Bl). 


C. OPTIMAL WEIGHTING OF MULTI-CHANNEL SPECTRA 

We use the 8 high-frequency differencing assemblies Q1 through W4 to generate the final combined spectrum. 
In this Appendix we show how we combine these spectra into a single estimate of the angular power spectrum. 
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The ultimate goal of the WMAP analysis is to produce a likelihood function for a set of cosmological parameters, 
{a}, given the data, Cj. Specifically 

V(a\q) = C(C i ,\Cf > (a))V(a}, (Cl) 

where V{a\C\) is the probability of {a} given the data, £(Cj|C| h (a)) is the likelihood of the data given the model, 
C^(a), and V(a) is the prior probability of the parameter set (Verde et al. 2003). To this end, we seek a combined 
spectrum, C/, that estimates the power spectrum in our sky, cf®, with the property that V{a\ C/) = V(a\ C]), and hence 

To estimate the combined spectrum, we approximate the likelihood function for the cross-power spectra as Gaus- 
sian 

-21n£(C*|C; h ) = £«3-C?wj) (Z^)% ( Cj , (C2) 

ij W 

where Cj is the spectrum with the best-fit source model subtracted, defined after equation (12), w) is the window 
function of spectrum i, defined after equation (6), and Efuii is the covariance matrix of the 28 cross-power spectra. 
Note that the treatment in the remainder of this section does not depend on any specific property of the covariance, 
so we use generic notation for readability. However, when we generate the WMAP first-year combined spectrum, the 
actual form of the covariance used at this step is (£f U n)Jj/, where E indicates the approximate covariance, obtained in 
§4, that has not yet had the effects of the foreground mask accounted for. 

We seek a spectrum Q such that 

-21n£(C,jCf ) = £(C, -Cf )&,. (Cj. -Cjf) = -21n£(C||C?), (C3) 

W 

where Qw is the inverse-covariance matrix of the combined spectrum which comes with the estimate of C/. Suppose, 
for simplicity, that E fu n is diagonal, (Efuii)//# = (E^n then it is straightforward to show that the deconvolved, 
weighted-average spectrum 

0 a — 1 

with 

Qi - 

U 

is the desired spectrum. Substituting equations (C4) and (C5) into equation (C3) produces equation (C2) up to a term 
which has a weak dependence on comological model, which we ignore. This combined spectrum is equivalent to the 
result we would obtain using the “optimal data compression” of Tegmark et al. (1997). 


(C4) 

(C5) 


If the inverse covariance matrix is not diagonal in /, it can be shown that the optimal combined spectrum is given 


by 


C/ = 




-£±ii 


E/-<2«'C ^ 


with 


U 


(C6) 

(C7) 


and where Cf d is a fiducial cosmological model which we take to be a flat ACDM model with = 0.021, £l c h 2 = 
0.129, ft tot = 1 , h = 0.68, n s = 1.2, and r = 0.2. While this model has parameters values that are substantially different 
than the best-fit WMAP model obtained (afterwards) by Spergel et al. (2003), the parameter degeneracies are such that 
Cf d is close to the best-fit model for / > 100 where this estimator is actually used (see §5.1). The combined spectrum 
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is optimal if the fiducial model chosen is the correct one; otherwise it is still unbiased but slightly sub-optimal (Gupta 
& Heavens 2002). 


D. CUT-SKY FISHER MATRIX 

The WMAP sky maps have nearly diagonal pixel-pixel noise covariance (Hinshaw et ai. 2003a) which greatly 
simplifies the properties of the power spectrum Fisher matrix. In this Appendix we present an analytic derivation of 
the effect of a sky cut and non-uniform pixel weighting. In D.l, we assume that we have an optimal estimator of 
the power spectrum. In the noise dominated limit, we can obtain an exact expression, while in the signal dominated 
limit, we need to approximate the signal correlation matrix to obtain an analytic expression. In D.2, we interpolate 
the Fisher matrix between the signal and noise-dominated regimes and show that it agrees with numerical estimates. 
In D.3, we estimate the power spectrum covariance matrix from Monte Carlo simulations of the sky and calibrate the 
interpolation formula. 


D.l. Cut-Sky Fisher Matrix: Analytic Evaluation 

We equate the Fisher matrix of the power spectrum to the curvature of the likelihood function, equation (A 16), 
then develop an approximate form for the covariance matrix, C = S + N, where S is the signal matrix and N is the noise 
matrix. We split the noise matrix into two pieces: a weight term and a mask term. In the limit that the pixel noise is 
diagonal, the weight term has the form 

(N w )ij = ~2^U = w $ij 0D1) 

where i and j are pixel indices, m is the number of observations of pixel i, ao is the rms noise of a single observation, 
and, by definition, w,- is the weight of pixel i. The mask, Af f -, is defined so that M\ equals 0 in pixels that are not used 
due to foreground contamination, and equals 1 otherwise. The noise matrix can then be written as the product of the 
two terms 

N" 1 = n; 3 m = mn; 1 = —Sij = w,$ y (D2) 

^0 

where hi - rii and wt = w, in the unmasked pixels, and are set to 0 otherwise. We thus define the covariance matrix over 
the full sky, which allows us to exploit the orthogonality properties of the spherical harmonics. Note that M 2 = M. 

The inverse of the full covariance matrix can now be written as 

C _1 = (S+N) -1 = N-^SN^+I)- 1 (D3) 

= M-^SM^+Nwr 1 . (D4) 

The covariance matrix has two limits. In the noise dominated limit, SN -1 <C I, 

C _1 — ► N" 1 . (D5) 

In the signal dominated limit, SM" 1 <C N w , only the mask alters the covariance matrix, so 

C' 1 -4 (D6) 

where we have set the inverse of the mask to zero where there is no data, i.e. M' 1 = M. 
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D.l.l. Cut-Sky Fisher Matrix in the Noise Dominated Limit 
In the noise dominated limit, the Fisher matrix, equation (A16), is approximately 


= ^tr(C _1 P;C' 1 P/<) -4 ^tr(N _1 P;N _1 P;') . 


Using 


(P/)/7 — op — 


DC (2/+1) 


~w { 


Pi(cosOij) = w, ^Y* m (pi)Yi m {pj) 


(D7) 


(D8) 


HJ ~ dC, 4tt 

m 

where 0^ is the angle between pixels i and j y Yi m (pi ) is a spherical harmonic evaluated in the direction of pixel i » and 
we have employed the addition theorem for spherical harmonics in the last equality. Then the Fisher matrix takes the 
form 

= X -w m . £ £ WMYUpj) WjY,t ml ( P j)Y rm , ( Pi ). (D9) 

~ mm 1 ij 

Now expand the weight array as 

Wi = W'tm Ylmipd (D 1 0) 




and substitute this into equation (D9) to yield 


F w = \ W W Y2 Y2 


Wi" m » Wi't’m'ft ^ Yl"m" (Pi)Y* m (pi)Ylm(p j ) (pj)Y Pm t (pj)Yr m > (pi). (Dll) 

mm' V’m"l’"m' n ij 


Since w,- was defined over the full sky, the sum over pixels may be expressed in terms of products of Wigner 3 -j 
symbols. As shown in Appendix E, we may use an orthogonality property of the 3 -j symbols to reduce the expression 
for the Fisher matrix to 


1 (21 + 1 )(2V + 1 ) ^ j (l V l" \ 

F "‘ 4,n;, X>(o o oj 


(D12) 


where W{ = |w/ m | 2 > and Q p is the solid angle per pixel. 


D. 1. 2. Cut- Sky Fisher Matrix in the Signal Dominated Limit 
In the signal dominated limit, the Fisher matrix, equation (A16), is approximately 

F? v = ^(CT’P/C-’P,.) -4 ] -tr (S-'P/S-’P,.) , 

where S~ l = MiSj-Mj. Using equation (D8) we can write F$, in the form 

F,f, = M J^j Y *JPj)YiM J2 Y Vm' (PpW'm' (Pd. 

” mm f i jk kp 

Substituting a normalized expression for in the pixel basis. 


(D13) 


(D14) 


(D15) 
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into equation (D1 4), one obtains a term 


' Cut Wl’t 

l" 1 1 m" i 


(D16) 


Since our mask cuts only ~ 15% of the sky, the coupling sum, M/7/w (pi)Yj!, m „ (pi), peaks very sharply at \V -/" | < 
/'. Therefore, one may approximate equation (D16) with 


ft 




(C/w/C/tW //) 1 / 2 

where, in the last equality, we have used the completeness relation for the spherical harmonics 

EC(p.)K/ m (/>,) = ft;%. 


/m 


With this approximation equation (D14) reduces to 


~ _ 
m/' - ~ 


i ft 2 


2 QC, 


;EE M Y* m {pi)Yi m (pj) Mj Y,: m , ( Pj )Y Vm , {pi). 


mm' ij 

Following the same steps outlined above for the derivation of equation (D12) yields 


ol 1 (2Z+1)(2Z'+1) 

t u t - ~ 


2 C\Cy 


47T 


?"'“(o 0 o) 


(D17) 


(D18) 


(D19) 


(D20) 


where Mi = Yl m \ m i™ i 2 > ancl spherical harmonic transform of the mask. Note that Mo = 4 tt/^ , where / S ky is 

the fraction of the sky that survives the mask. 


D.2. Interpolating the Cut-Sky Fisher Matrix 


We can combine the two limiting cases of the Fisher matrix, obtained in the previous section, into a single 
expression 


17 mask _ 
t lV “ 


V(2Z+1)(2Z'+1)/ 2 
2 (Ci+N,)(C,.+N l .) 


(D21) 


where TV/ = n p cro/(«^uv) is the deconvolved noise power spectrum. Using equations (D12) and (D20) derived in the 
previous section, Fw can be expressed in the low l (signal-dominated) and high / (noise-dominated) limits as 


and 


rS _ 
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V(2/+l)(2/' + l) 


4tt/: 2 
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(D23) 


By construction, these matrices are normalized to 1 on the diagonal, F ; ? = F* = 1, since Mq = 47r/ 2 kv and VV 0 = 

47r /4y«L/ cr 0- 
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We have computed directly, for selected / values, by evaluating Yi m C~ x Yy m > using the pre-conditioner code 
described by Oh et al. (1999). We find that deviations between the numerical and analytical results are consistent with 
numerical noise in the Fisher matrix estimate. Figure 14 shows the estimated Fisher matrix for the Kp2 cut in the noise 
and the signal dominated limits. 


We interpolate between the two regimes with the following expression: 

.mask 

F “' (G+NiXCe+Nr) 


(D24) 


Note that the form of the Fisher matrix primarily depends on A/, but is weakly dependent on / - there is more coupling 
in the noise dominated limit. 


D.3. Power Spectrum Covariance: Monte Carlo Evaluation 

As discussed in Appendix A. 1 .2 we compute the Q using three different pixel weightings. The uniform weighting 
and the inverse-noise inverse weighting are optimal in the signal-dominated regime and the noise-dominated regime, 
respectively. In between these limits the transitional weighting performs better. In order to determine which ranges in / 
correspond to which regimes, and to calibrate our ansatz for the covariance matrix, equation (D21), we proceed as fol- 
lows. Using 100,000 Monte Carlo simulations of signal plus noise (with WMAPnoise levels and symmetrized beams), 
we compute the diagonal elements of the covariance matrix for the three different weighting schemes, evaluated with 
the Kp2 sky cut. Denote these estimates Df m 

We find that the uniform weighting produces the smallest Df m below / = 200. Inverse-noise weighting is the best 
above / = 450, and the transitional weighting produces the lowest variance in between. In each of these regimes we 
use the resulting Df m to calibrate our ansatz for the diagonal elements of the covariance matrix 

Dl ‘WTi^ c '* Nf)2 ’ <D25) 

as illustrated in Verde et al. (2003). These calibrations produce a smooth correction to equation (D24) of at most 6%. 
No correction at all is required in the signal dominated regime. 


E. SOME USEFUL PROPERTIES OF SPHERICAL HARMONICS 

The derivation of the form of the Fisher matrix in the signal and noise dominated limits led to expressions which 
included a term 

53 E ( PiKm(Pi)Ylm(Pj ) (Pj) ( Pj)Y Vm . ( Pi ), (El ) 

mm 1 i j 

where the sum is effectively a double integral over the full sky. This can be evaluated in terms of the Wigner 3 -j 
symbols, defined as 

t\m' . /(2J+ 1)(2/'+ l)(2/" + l) ( l V l"\(l V l" \ 
j dn„Y Im (n)Y l , ml (n)Y l , lm „(n) = (-1) y — ( 0 0 oJ(m -m' m")' (E2) 

Substituting this into equation (El) gives 


4-7T 
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We simplify this using the Wigner 3- j orthogonality condition 


(E3) 



$m n m ,n 5(1, l f , l”), 


(E4) 


where 8(1, V,V ') = 1 for \l~V\ < V < l+l f and is 0 otherwise. This reduces equation (El) to 

(27+l)(2/'+l) /' I l ' l" \ 2 

4nQp ( 0 0 0 / ’ (E5) 

where the factor of Q 2 p accounts for the fact that equation (El ) is a double sum over pixels, instead of a double integral. 
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Fig. 1. — The effective noise as a function of / for the 8 differencing assemblies used in the combined power spectrum 
analysis. These spectra were computed from end-end simulations of noise maps, as discussed in §2.2. For illustration, 
the spectra shown here were computed using the quadratic estimator with uniform pixel weights. The actual noise 
model uses three separate weighting schemes in three separate l ranges, and thus has discontinuities where the effective 
noise level changes. The weights are defined in Appendix A. 1.2. 
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Fig. 2. — The full set of individual cross-power spectra for l < 500, computed from the 8 high frequency differencing 
assemblies Q1 through W4, 28 spectra in ail. The spectra were evaluated from the uncorrected sky maps (no Galaxy 
model subtracted) using the Kp2 sky cut with uniform weighting. The data are plotted in color by effective frequency 
y/uiVj with red corresponding to 41 GHz and blue to 94 GHz. The top panel shows a very robust measurement of the 
first acoustic peak with a maximum near / ~ 220. There is also a clear indication of the rise to a second peak at / ~ 540 
as discussed in §6. The bottom panel shows the ratio of each channel to the combined spectrum presented in §5. This 
clearly shows the residual foreground emission due to diffuse Galactic radio emission at low / and to point sources at 
higher /. The level of contamination, which is strongest at Q band, is consistent with the expected level of foreground 
emission. See Figure 3 for the spectra after foreground subtraction. 
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Fig. 3. — The same set of cross-power spectra as shown in Figure 2. Here, the foreground model discussed in §3.1 
has been subtracted from each channel. The bottom panel shows the ratio of each of the 28 cross-power spectra to 
the combined spectrum presented in §5. Aside from a ~10% discrepancy in the Q band data at / < 20, the data 
are consistent with each other to the sensitivity limits of the individual spectra. Because the WMAP data are not 
sensitivity limited at low /, we use only V and W band data in the final combined spectrum for / < 100 (see §5) to 
minimize residual Galactic contamination. 
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Fig. 4. — The same set of cross-power spectra as shown in Figure 3, but showing the low / spectrum unbinned. The 
agreement between the individual spectra is striking. The low value of the quadrupole moment, C 2 , that was first seen 
by COBE-DMR is also seen in the WMAPdata. The steep, nearly linear rise in the spectrum from / = 2 to 5 translates 
to a near absence of power in the angular correlation function at separations larger than ~60° (Spergel et al. 2003; 
Bennett et al. 2003b). This was also seen in the COBE-DMR data, but it is now clear that this is not due to Galaxy 
modeling errors. The bottom panel shows the fractional rms among the 28 WMAP cross-power spectra in black, while 
the red curve shows the same statistic averaged over an ensemble of 1000 Monte Carlo realizations. Based on this we 
estimate the measurement error on the combined spectrum to be <2-5% for / < 100. 
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Fig. 5. — The auto-power spectrum of the combined Q+V+W map evaluated with the three weighting schemes defined 
in Appendix A. 1 .2. In each case, the spectrum was computed over the entire / range, and black shows uniform weights, 
red shows inverse-noise weights, and green shows transitional weights. The agreement is excellent. 
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Fig. 6. — This figure compares the auto-power spectrum computed from the combined Q+V+W map (grey) to the 
optimally combined cross-power spectrum (black). The close agreement indicates that the noise properties of the 
first-year WMAP data are well understood 
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Fig. 7. — A comparison of the power spectrum computed with data from the ecliptic plane (black) vs. data from the 
ecliptic poles (grey). Note that some of the “bite” features that appear in the combined spectrum are not robust to data 
excision. There is also no evidence that beam ellipticity, which would be more manifest in the plane than in the poles, 
systematically biases the spectrum. This is consistent with estimates of the effect given by Page et al. (2003a). 




-36- 



Angular scale (deg) 

90 2 0.5 0.2 



10 40 100 200 400 800 1 400 

Multipole moment L 


Fig. 8. — The final angular power spectrum, /(/ + l)C//27r, obtained from the 28 cross-power spectra, as described in 
§5. The data are plotted with la measurement errors only which reflect the combined uncertainty due to noise, beam, 
calibration, and source subtraction uncertainties. The solid line shows the best-fit ACDM model from Spergel et al. 
(2003). The grey band around the model is the la uncertainty due to cosmic variance on the cut sky. For this plot, both 
the model and the error band have been binned with the same boundaries as the data, but they have been plotted as a 
splined curve to guide the eye. On the scale of this plot the unbinned model curve would be virtually indistinguishable 
from the binned curve except in the vicinity of the third peak. 
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Fig. 9. — A compilation of recent CMB power spectrum measurements compared to the best-fit ACDM model from the 
first-year WMAPdata. The data points include noise and cosmic variance uncertainty (but not calibration uncertainty) 
thus we omit the cosmic variance band from the model curve in the Figure. On average, the pre- WMAPdata agree 
well with the WMAP power spectrum. 
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Fig. 10. — The WMAP combined power spectrum, in black, compared to a compilation of all CMB data published 
prior to WMAP from Wang et al. (2002a), in red. The WMAP data are plotted with cosmic variance plus measurement 
uncertainties here in order to facilitate a comparison with the compiled data which is reported in this way. The data 
agree well on COBE scales, / < 20, (but note that the WMAP cosmic variance errors are computed from the best-fit 
model rather than the data, thus they appear larger than the COBE errors at the quadrupole). However, the overall 
normalization of the WMAP spectrum is ~10% higher on smaller scales. 
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Fig. 1 1 . — The WMAP combined power spectrum compared to the locus of predicted spectra, in red, based on a joint 
analysis of pre-WMAPCMB data and 2dFGRS large-scale structure data (Percival et al. 2002). As in Figure 8, the 
WMAP data are plotted with measurement uncertainties, and the best-fit ACDM model (Spergel et al. 2003) is plotted 
with ala cosmic variance error band. The locus of predicted spectra lie systematically below the WMAP data at 
intermediate /. 
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Fig. 12. — A comparison of two power spectra computed from a single V band map. The black points result from 
the maximum likelihood method, the red points from the quadratic estimator computed with uniform pixel weighting. 
In both cases a noise model has been assumed to treat the noise bias. At higher l, the two spectra would be nearly 
identical since both impose inverse noise weighting on the data. 
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Fig. 13. — The quadratic estimator discussed in §A.l has been extensively tested with Monte Carlo simulations of 
first-year WMAP data. The top panel shows a model spectrum in red, and the mean of 600 realizations of the cross- 
power spectrum, computed from the V and W band-averaged maps, in black. The red error bars show' the uncertainty 
in the mean in bins with A l = 20. The cross-power spectrum estimator is unbiased. The bottom panel is similar, but 
with auto-power spectra computed from the W band-averaged map. Here the noise bias term was estimated directly 
from the high / tail of the computed spectrum. This simple auto-power spectrum estimator is unbiased if the noise is 
white, and is suitable for V band data, but see Figure 1 for examples of W band noise spectra. 
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Fig. 14. — Slices of the Fisher (or curvature) matrix normalized as Fw / VFiFrr - plotted vs. A l = 1-1', from the Fisher 
matrix for the combined spectrum. Black is / = 30, red is / = 300, green is / = 600. 
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Table 1 . WMAP Power Spectrum Noise Model 3 


Parameter Q1 Q2 VI V2 W1 W2 W3 W4 


/ < 200 


co(xl0 2 ) 

-6.604 

-6.971 

-6.823 

-6.674 

-7.174 

-7.439 

-7.044 

-7.989 

cj (x 10' 1 ) 

1.934 

4.244 

2.497 

1.395 

3.060 

3.721 

1.459 

4.820 

c 2 (x]0 4 ) 

-3.798 

-8.521 

-5.203 

-1.483 

-6.021 

-6.901 

-2.677 

-5.870 

C3(Xl0 5 ) 

2.549 

5.921 

3.715 

0.242 

4.147 

4.540 

1.750 

2.167 


200 </ < 450 


co(xl0 2 ) 

-6.143 

-6.268 

-6.320 

-6.192 

-6.618 

-6.833 

-6.768 

-6.849 

O(xl0 4 ) 

-1.334 

1.341 

-0.800 

-1.351 

0.769 

2.582 

0.885 

2.443 

/ > 450 

c 0 (xl0 2 ) 

-6.202 

-6.167 

-6.338 

-6.264 

-6.564 

-6.704 

-6.688 

-6.772 

ci(xl0 5 ) 

3.288 

4.085 

3.327 

5.979 

7.342 

12.970 

8.851 

24.430 


a Best-fi t coeffi cients for the noise model in equation (8). The units of the output noise are ml£ , 
thermodynamic temperature. 










